
!<    @file     This file contains the LBFGS algorithm and supporting routines
!>
!>@brief        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
!>                          JORGE NOCEDAL
!>                        *** July 1990 ***
!>
!> 
!>     This subroutine solves the unconstrained minimization problem
!> 
!>                      min F(x),    x= (x1,x2,...,xN),
!>
!>      using the limited memory BFGS method. The routine is especially
!>      effective on problems involving a large number of variables. In
!>      a typical iteration of this method an approximation Hk to the
!>      inverse of the Hessian is obtained by applying M BFGS updates to
!>      a diagonal matrix Hk0, using information from the previous M steps.
!>      The user specifies the number M, which determines the amount of
!>      storage required by the routine. The user may also provide the
!>      diagonal matrices Hk0 if not satisfied with the default choice.
!>     The algorithm is described in "On the limited memory BFGS method
!>      for large scale optimization", by D. Liu and J. Nocedal,
!>      Mathematical Programming B 45 (1989) 503-528.
!> 
!>      The user is required to calculate the function value F and its
!>      gradient G. In order to allow the user complete control over
!>      these computations, reverse  communication is used. The routine
!>      must be called repeatedly under the control of the parameter
!>      IFLAG. 
!>
!>      The steplength is determined at each iteration by means of the
!>      line search routine MCVSRCH, which is a slight modification of
!>      the routine CSRCH written by More' and Thuente.
!>
!>    ON THE DRIVER: The program that calls LBFGS must contain the declaration:
!>
!>                       EXTERNAL LB2
!>
!>    LB2 is a BLOCK DATA that defines the default values of several
!>    parameters described in the COMMON section. 


!>    COMMON:
!> 
!>     The subroutine contains one common area, which the user may wish to
!>   reference:
!> 


!>  MACHINE DEPENDENCIES
!>
!>        The only variables that are machine-dependent are XTOL,
!>        STPMIN and STPMAX.
!> 
!>
!>  GENERAL INFORMATION
!> 
!>    Other routines called directly:  DAXPY, DDOT, LB1, MCSRCH
!> 
!>    Input/Output  :  No input; diagnostic messages on unit MP and
!>                     error messages on unit LP.
!> 
!> 


!> @param    N       is an INTEGER variable that must be set by the user to the number of variables. It is not altered by the routine. Restriction: N>0.

!> @param    M       is an INTEGER variable that must be set by the user to the number of corrections used in the BFGS update. It is not altered by the routine. Values of M less than 3 are not recommended; large values of M will result in excessive computing time. 3<= M <=7 is recommended. Restriction: M>0.

!> @param    ITER       numero da iteracao (acrescentado como parâmetro por franciane em 8/7/2011)

!> @param    X       is a DOUBLE PRECISION array of length N. On initial entry it must be set by the user to the values of the initial estimate of the solution vector. On exit with IFLAG=0, it contains the values of the variables at the best point found (usually a solution).

!> @param    F       is a DOUBLE PRECISION variable. Before initial entry and on a re-entry with IFLAG=1, it must be set by the user to contain the value of the function F at the point X.

!> @param    G       is a DOUBLE PRECISION array of length N. Before initial entry and on a re-entry with IFLAG=1, it must be set by the user to contain the components of the gradient G at the point X.

!> @param    DIAGCO  is a LOGICAL variable that must be set to .TRUE. if the user  wishes to provide the diagonal matrix Hk0 at each iteration. Otherwise it should be set to .FALSE., in which case  LBFGS will use a default value described below. If DIAGCO is set to .TRUE. the routine will return at each iteration of the algorithm with IFLAG=2, and the diagonal matrix Hk0  must be provided in the array DIAG.


!> @param    DIAG    is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE., then on initial entry or on re-entry with IFLAG=2, DIAG it must be set by the user to contain the values of the diagonal matrix Hk0.  Restriction: all elements of DIAG must be positive.

!> @param    IPRINT  is an INTEGER array of length two which must be set by the user.
!>            IPRINT(1) specifies the frequency of the output:
!>               IPRINT(1) < 0 : no output is generated,
!>               IPRINT(1) = 0 : output only at first and last iteration,
!>               IPRINT(1) > 0 : output every IPRINT(1) iterations.
!>            IPRINT(2) specifies the type of output generated:
!>               IPRINT(2) = 0 : iteration count, number of function evaluations, function value, norm of the gradient, and steplength,
!>               IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of  variables and  gradient vector at the initial point,
!>               IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of variables,
!>               IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.


!> @param    EPS     is a positive DOUBLE PRECISION variable that must be set by the user, and determines the accuracy with which the solution is to be found. The subroutine terminates when ||G|| < EPS max(1,||X||), where ||.|| denotes the Euclidean norm.

!> @param    XTOL    is a  positive DOUBLE PRECISION variable that must be set by the user to an estimate of the machine precision (e.g. 10**(-16) on a SUN station 3/60). The line search routine will terminate if the relative width of the interval of uncertainty is less than XTOL.

!> @param    W       is a DOUBLE PRECISION array of length N(2M+1)+2M used as workspace for LBFGS. This array must not be altered by the user.

!> @param    IFLAG   is an INTEGER variable that must be set to 0 on initial entry to the subroutine. 
!> A return with IFLAG<0 indicates an error, and IFLAG=0 indicates that the routine has terminated without detecting errors. On a return with IFLAG=1, the user must evaluate the function F and gradient G. On a return with IFLAG=2, the user must provide the diagonal matrix Hk0.
!>            The following negative values of IFLAG, detecting an error, are possible:
!>             IFLAG=-1  The line search routine MCSRCH failed. The parameter INFO provides more detailed information (see also the documentation of MCSRCH):
!>                     INFO = 0  IMPROPER INPUT PARAMETERS. 
!>                     INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL.
!>                      INFO = 3  MORE THAN 20 FUNCTION EVALUATIONS WERE REQUIRED AT THE PRESENT ITERATION.
!>                      INFO = 4  THE STEP IS TOO SMALL.
!>                      INFO = 5  THE STEP IS TOO LARGE.
!>                      INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS. THERE MAY NOT BE A STEP WHICH SATISFIES THE SUFFICIENT DECREASE AND CURVATURE CONDITIONS. TOLERANCES MAY BE TOO SMALL.
!>             IFLAG=-2  The i-th diagonal element of the diagonal inverse Hessian approximation, given in DIAG, is not positive.
!>             IFLAG=-3  Improper input parameters for LBFGS (N or M are not positive).
!> 
!> @param      MP  is an INTEGER variable with default value 6. It is used as the
!>        unit number for the printing of the monitoring information
!>        controlled by IPRINT.
!> 
!> @param      LP  is an INTEGER variable with default value 6. It is used as the
!>        unit number for the printing of error messages. This printing
!>        may be suppressed by setting LP to a non-positive value.
!> 
!> @param     GTOL is a DOUBLE PRECISION variable with default value 0.9, which
!>        controls the accuracy of the line search routine MCSRCH. If the
!>        function and gradient evaluations are inexpensive with respect
!>        to the cost of the iteration (which is sometimes the case when
!>        solving very large problems) it may be advantageous to set GTOL
!>        to a small value. A typical small value is 0.1.  Restriction:
!>        GTOL should be greater than 1.D-04.
!> 
!>  @param     STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
!>        specify lower and uper bounds for the step in the line search.
!>        Their default values are 1.D-20 and 1.D+20, respectively. These
!>        values need not be modified unless the exponents are too large
!>        for the machine being used, or unless the problem is extremely
!>       badly scaled (in which case the exponents should be increased).
!> 
C
C     THE WORK VECTOR W IS DIVIDED AS FOLLOWS:
C     ---------------------------------------
C     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND
C         OTHER TEMPORARY INFORMATION.
C     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.
C     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED
C         IN THE FORMULA THAT COMPUTES H*G.
C     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH
C         STEPS.
C     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M
C         GRADIENT DIFFERENCES.
C
C     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A
C     CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT.
C


      SUBROUTINE lmLAMBDALBFGS(N,M,ITER,X,F,F1,F2,G,G1,G2,lambda,
     &      DIAGCO,DIAG,IPRINT,EPS,XTOL,Q,Q1,Q2,RHO,ALPHA,S,
     &      Y,Y1,Y2,X0,IFLAG)
C
      implicit none
      INTEGER N
      INTEGER M
      INTEGER IPRINT(2)
      INTEGER IFLAG
      DOUBLE PRECISION X(N)
      DOUBLE PRECISION G(N)
      DOUBLE PRECISION G1(N)
      DOUBLE PRECISION G2(N)
      DOUBLE PRECISION lambda
      DOUBLE PRECISION DIAG(N)

!       DOUBLE PRECISION W(N*(2*M+1)+2*M)
      DOUBLE PRECISION Q(N)
      DOUBLE PRECISION Q1(N)
      DOUBLE PRECISION Q2(N)
      DOUBLE PRECISION RHO(M)
      DOUBLE PRECISION ALPHA(M)
      DOUBLE PRECISION S(N*M)!< PASSOS DAS ÚLTIMAS M ITERAÇÕES
      DOUBLE PRECISION Y(N*M)!< DIFERENCAS ENTRE GRADIENTES DAS ÚLTIMAS  ITERAÇÕES
      DOUBLE PRECISION Y1(N*M)!< DIFERENCAS ENTRE GRADIENTES DAS ÚLTIMAS  ITERAÇÕES, parcela do resíduo 
      DOUBLE PRECISION Y2(N*M)!< DIFERENCAS ENTRE GRADIENTES DAS ÚLTIMAS  ITERAÇÕES , parcela do termo de regularização    
      DOUBLE PRECISION X0(N)!< X da iteração anterior  


      DOUBLE PRECISION F
      DOUBLE PRECISION F1
      DOUBLE PRECISION F2
      DOUBLE PRECISION F1min
      DOUBLE PRECISION EPS
      DOUBLE PRECISION XTOL
      LOGICAL DIAGCO

C
C 
C 

         COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX

      DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN,
     .                 STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM
      INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO,
     .        BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN
      LOGICAL FINISH
C
      SAVE
      DATA ONE,ZERO/1.0D+0,0.0D+0/


C
C     INITIALIZE
C     ----------
C

      IF(IFLAG.EQ.0) GO TO 10
      GO TO (172,100) IFLAG
  10  ITER= 0

      F1min=F1



      IF(N.LE.0.OR.M.LE.0) GO TO 196
      IF(GTOL.LE.1.D-04) THEN
	  IF(LP.GT.0) WRITE(LP,245)
	  GTOL=9.D-01
      ENDIF
      NFUN= 1
      POINT= 0
      FINISH= .FALSE.
      IF(DIAGCO) THEN
	  DO I=1,N
	      IF (DIAG(I).LE.ZERO) GO TO 195
	  ENDDO
      ELSE
	  DO I=1,N
	      DIAG(I)= 1.0D0
	  ENDDO
      ENDIF

      ISPT= N+2*M
      IYPT= ISPT+N*M     
      DO I=1,N
	  !  50   W(ISPT+I)= -G(I)*DIAG(I)
	  S(I)= -G(I)*DIAG(I)
      ENDDO
      GNORM= DSQRT(DDOT(N,G,1,G,1))
      STP1= ONE/GNORM
C
C     PARAMETERS FOR LINE SEARCH ROUTINE
C     
      FTOL= 1.0D-4!< (original,ftol=1.0D-4)
      MAXFEV= 20 !< número maximo de avaliações por iteração (original, maxfev=20)
C
      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
     *                     GNORM,N,M,X,F,G,STP,FINISH)
C
C    --------------------
C     MAIN ITERATION LOOP
C    --------------------
C
 80   ITER= ITER+1

      do i=1,N
	  X0(i) = X(i)
      enddo

      INFO=0
      BOUND=ITER-1
      IF(ITER.EQ.1) GO TO 165
      IF (ITER .GT. M)BOUND=M
C
!          YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1)
         YS= DDOT(N,Y(NPT+1),1,S(NPT+1),1)
      IF(.NOT.DIAGCO) THEN
!          YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1)
	  YY= DDOT(N,Y(NPT+1),1,Y(NPT+1),1)
	  DO I=1,N
	      DIAG(I)= YS/YY
	  ENDDO
      ELSE
	  IFLAG=2
	  RETURN
      ENDIF
 100  CONTINUE
      IF(DIAGCO) THEN
	  DO I=1,N
	      IF (DIAG(I).LE.ZERO) GO TO 195
	  ENDDO
      ENDIF
C
C     COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
C     "Updating quasi-Newton matrices with limited storage",
C     Mathematics of Computation, Vol.24, No.151, pp. 773-782.
C     ---------------------------------------------------------
C
      CP= POINT
      IF (POINT.EQ.0) CP=M
!       W(N+CP)= ONE/YS
      rho(CP)= ONE/YS
      DO I=1,N
	  !  112  W(I)= -G(I)
	  Q(I)= -G(I)
      ENDDO
      CP= POINT

      DO I= 1,BOUND
	  CP=CP-1
	  IF (CP.EQ. -1)CP=M-1

	  !SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1)
	  SQ= DDOT(N,S(CP*N+1),1,Q,1)

	  INMC=N+M+CP+1
	  IYCN=IYPT+CP*N

	  !W(INMC)= W(N+CP+1)*SQ
	  ALPHA(CP+1)= RHO(CP+1)*SQ

	  !CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1)
	  CALL DAXPY(N,-ALPHA(CP+1),Y(CP*N+1),1,Q,1)
      ENDDO
C
      DO I=1,N
	  ! 130  W(I)=DIAG(I)*W(I)
	  Q(I)=DIAG(I)*Q(I)
      ENDDO
C
      DO I=1,BOUND
	  !YR= DDOT(N,W(IYPT+CP*N+1),1,W,1)
	  YR= DDOT(N,Y(CP*N+1),1,Q,1)

	  !BETA= W(N+CP+1)*YR
	  BETA= RHO(CP+1)*YR

	  INMC=N+M+CP+1

	  !BETA= W(INMC)-BETA
	  BETA= ALPHA(CP+1)-BETA

	  ISCN=ISPT+CP*N

	  !CALL DAXPY(N,BETA,W(ISCN+1),1,W,1)
	  CALL DAXPY(N,BETA,S(CP*N+1),1,Q,1)

	  CP=CP+1
	  IF (CP.EQ.M)CP=0
      ENDDO
C
C     STORE THE NEW SEARCH DIRECTION
C     ------------------------------
C
      DO I=1,N
	  !160   W(ISPT+POINT*N+I)= W(I)
	  S(POINT*N+I)= Q(I)
      ENDDO
C
C     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION 
C     BY USING THE LINE SEARCH ROUTINE MCSRCH
C     ----------------------------------------------------
 165  NFEV=0
      STP=ONE
      IF (ITER.EQ.1) STP=STP1
      DO I=1,N
	  !  170  W(I)=G(I)
	  Q(I)=G(I)
	  Q1(I)=G1(I)
	  Q2(I)=G2(I)
      ENDDO
 172  CONTINUE
!       CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,
!      *            XTOL,MAXFEV,INFO,NFEV,DIAG)
      CALL MCSRCH(N,X,F,G,S(POINT*N+1),STP,FTOL,
     *            XTOL,MAXFEV,INFO,NFEV,DIAG)
      IF (INFO .EQ. -1) THEN
	  IFLAG=1
	  RETURN
      ENDIF
      IF (INFO .NE. 1) GO TO 190


! ! atualização do parâmetro da regularização:

!       if(F1.lt.F1min.and. F2.ne.0.0d0)then
! 	  F1min = F1
! 	  lambda = F1/F2
!       endif

      if(F1.lt.lambda*F2)then!equivalente a F1/F2.lt.lambda, mas na possibilidade de F2 ser zero, escrevi deste jeito
	  lambda = F1/F2
      endif
      write(*,*) 'LAMBDALBFGS    lambda=',lambda



      NFUN= NFUN + NFEV
C
C     COMPUTE THE NEW STEP AND GRADIENT CHANGE 
C     -----------------------------------------
C
      NPT=POINT*N
      DO I=1,N
	  ! W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
	  S(NPT+I)= STP*S(NPT+I)

	  !  175  W(IYPT+NPT+I)= G(I)-W(I)
! 	  Y(NPT+I)= G(I)-Q(I)
	  Y1(NPT+I)= G1(I)-Q1(I) ! isso só é possivel pq lambda é o lambdaatual
	  Y2(NPT+I)= G2(I)-Q2(I)

      ENDDO

      do i=1,n*m
	  Y(I)= Y1(I) + lambda* Y2(I) ! mas daqui pra frente já usa o lambda novo para todos os Y anteriores
      enddo

      POINT=POINT+1
      IF (POINT.EQ.M)POINT=0
C
C     TERMINATION TEST
C     ----------------
C
      GNORM= DSQRT(DDOT(N,G,1,G,1))
      XNORM= DSQRT(DDOT(N,X,1,X,1))
      XNORM= DMAX1(1.0D0,XNORM)
      IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
C
      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
     *               GNORM,N,M,X,F,G,STP,FINISH)
      IF (FINISH) THEN
         IFLAG=0
         RETURN
      ENDIF
      GO TO 80
C
C     ------------------------------------------------------------
C     END OF MAIN ITERATION LOOP. ERROR EXITS.
C     ------------------------------------------------------------
C
 190  IFLAG=-1
      IF(LP.GT.0) WRITE(LP,200) INFO
      RETURN
 195  IFLAG=-2
      IF(LP.GT.0) WRITE(LP,235) I
      RETURN
 196  IFLAG= -3
      IF(LP.GT.0) WRITE(LP,240)
C
C     FORMATS
C     -------
C
 200  FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE'
     .          ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN'
     .          ' OF LINE SEARCH: INFO= ',I2,/
     .          ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,
     .          ' OR INCORRECT TOLERANCES')
 235  FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,
     .       ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
 240  FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',
     .       ' ARE NOT POSITIVE)')
 245  FORMAT(/'  GTOL IS LESS THAN OR EQUAL TO 1.D-04',
     .       / ' IT HAS BEEN RESET TO 9.D-01')
      RETURN
      END
C
C     LAST LINE OF SUBROUTINE LBFGS
C
C